##### Create all tables for SI

library(purrr)
library(texreg)
library(xtable)
library(readxl)
library(openintro)
library(rlang)
library(broom)
library(foreign)
library(purrr)
library(plm)
library(texreg)
library(foreign)
library(readxl)
library(openintro)
library(lubridate)
library(tidyverse)

source("lib/model-functions-HPBM.R") #these have all the different models we need to run
source("lib/fix-names.R")

## Create Supplementary Table 1 (Main BG Results)
load('data/results-BG-for-tables.RData')

table2.A <- list(BG_OLS_results[[3]], BG_OLS_results[[2]],BG_OLS_results[[1]],
                 first_stage_results[[1]], first_stage_results[[3]], first_stage_results[[2]],
                 BG_results[[1]], BG_county_results[[1]], BG_agency_results[[1]])

table2.B <- list(BG_results[[2]], BG_county_results[[2]], BG_agency_results[[2]],
                 BG_results[[3]], BG_county_results[[3]], BG_agency_results[[3]],
                 BG_results[[4]], BG_county_results[[4]], BG_agency_results[[4]])

table2.C <- list(BG_results[[5]], BG_county_results[[5]], BG_agency_results[[5]],
                 BG_results[[6]], BG_county_results[[6]], BG_agency_results[[6]],
                 BG_results[[7]], BG_county_results[[7]], BG_agency_results[[7]])

coef.map <- list(milex_iv = "Military exp. IV",
                 ltotal_cost = "Lagged Total Aid",
                 "log(1 + sumallLESOaidBG_lag)" = "Lagged Total Aid",
                 "`log(1 + sumallLESOaidBG_lag)(fit)`" = "Lagged Total Aid",
                 "`ltotal_cost(fit)`" = "Lagged Total Aid")

texreg(table2.A,
       file = "tables/BGtable2-partA-new.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.B,
       file = "tables/BGtable2-partB-new.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.C,
       file = "tables/BGtable2-partC-new.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Tables 2-4 (BG First Stage and Reduced Form)
load('data/results-BG-for-tables.RData')

bg <- list(first_stage_results$BG, BG_red_results[[1]], BG_red_results[[2]], BG_red_results[[3]],
           BG_red_results[[4]], BG_red_results[[5]], BG_red_results[[6]], BG_red_results[[7]])

county <- list(first_stage_results$rep_county, BG_county_red_results[[1]], BG_county_red_results[[2]], BG_county_red_results[[3]],
           BG_county_red_results[[4]], BG_county_red_results[[5]], BG_county_red_results[[6]], BG_county_red_results[[7]])

agency <- list(first_stage_results$rep_agency, BG_agency_red_results[[1]], BG_agency_red_results[[2]], BG_agency_red_results[[3]],
               BG_agency_red_results[[4]], BG_agency_red_results[[5]], BG_agency_red_results[[6]], BG_agency_red_results[[7]])

r.model.names <- c("First Stage", "Total Crime Rate", "Homicide", "Robbery", "Assault", 
                   "Burglary", "Larceny", "Vehicle Theft")

r.coef.names <- c("Percent Poverty", "log(Median Income)", "Unemployment Rate",
                  "log(Population)","Share Males", "Share Blacks", "Share Age 15-19",
                  "Share Age 20-24", "Share Age 25-29", "Share Age 30-34", "BG IV")

texreg(bg,
       single.row = FALSE,
       custom.model.names = r.model.names,
       custom.coef.names = r.coef.names,
       reorder.coef = c(11,seq(1,10)),
       caption = "First Stage and Reduced Form estimates of BG Results",
       label = "tab:reducedBG",
       script.size=TRUE,
       file="tables/BG-BG-reduced-form.tex",
       ci.force=T,
       ci.test = NULL)

texreg(county,
       single.row = FALSE,
       custom.model.names = r.model.names,
       custom.coef.names = r.coef.names[-10],
       reorder.coef = c(10,seq(1,9)),
       caption = "First Stage and Reduced Form estimates of County Level Replication Results",
       label = "tab:reducedCounty",
       script.size=TRUE,
       file="tables/BG-county-reduced-form.tex",
       ci.force=T,
       ci.test = NULL)

texreg(agency,
       single.row = FALSE,
       custom.model.names = r.model.names,
       custom.coef.names = r.coef.names[-10],
       reorder.coef = c(10,seq(1,9)),
       caption = "First Stage and Reduced Form estimates of Agency Level Replication Results",
       label = "tab:reducedAgency",
       script.size=TRUE,
       file="tables/BG-agency-reduced-form.tex",
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 5 (Main HPBM Results)
load('data/results-HPBM-for-tables.RData')

partA.1 <- list(HPBM_items_results[[1]], HPBM_county_items_results[[1]], HPBM_agency_items_results[[1]],
                HPBM_items_results[[2]], HPBM_county_items_results[[2]], HPBM_agency_items_results[[2]])

partA.2 <- list(HPBM_value_results[[1]], HPBM_county_value_results[[1]], HPBM_agency_value_results[[1]],
                HPBM_value_results[[2]], HPBM_county_value_results[[2]], HPBM_agency_value_results[[2]])

partB.1 <- list(HPBM_items_results[[3]], HPBM_county_items_results[[3]], HPBM_agency_items_results[[3]],
                HPBM_items_results[[4]], HPBM_county_items_results[[4]], HPBM_agency_items_results[[4]])

partB.2 <- list(HPBM_value_results[[3]], HPBM_county_value_results[[3]], HPBM_agency_value_results[[3]],
                HPBM_value_results[[4]], HPBM_county_value_results[[4]], HPBM_agency_value_results[[4]])

coef.map1 <- list("`log_sum_items_lag(fit)`" = "log items$_{t-1}$",
                  "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "log items$_{t-1}$")
coef.map2 <- list("`log_sum_value_ttl_lag(fit)`" = "log value$_{t-1}$",
                  "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "log value$_{t-1}$")

texreg(partA.1,
       file = "tables/HPBMtable8-partA1-new.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partA.2,
       file = "tables/HPBMtable8-partA2-new.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.1,
       file = "tables/HPBMtable8-partB1-new.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.2,
       file = "tables/HPBMtable8-partB2-new.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Tables 6-8 (HPBM First Stage and Reduced Form)
load('data/results-HPBM-for-tables.RData')

harris.models <- list(HPBM_first_stage_results[[5]], HPBM_red_items_results[[1]], HPBM_red_items_results[[2]],
                      HPBM_red_items_results[[3]], HPBM_red_items_results[[4]],
                      HPBM_first_stage_results[[6]], HPBM_red_value_results[[1]], HPBM_red_value_results[[2]],
                      HPBM_red_value_results[[3]], HPBM_red_value_results[[4]])

h.model.names <- c("First Stage", "Homicide", "Robbery", "Assault", "Vehicle Theft",
                   "First Stage", "Homicide", "Robbery", "Assault", "Vehicle Theft")

h.coef.names <- c("Percent Poverty","Unemployment Rate","log(Population)", 'Econmiss',
                  'Murder Arrests$_{t-1}$','Rape Arrests$_{t-1}$','Robbery Arrests$_{t-1}$',
                  'Agv. Assault Arrests$_{t-1}$', 'Burglary Arrests$_{t-1}$', 'Assault Arrests$_{t-1}$','Weapons Arrests$_{t-1}$',
                  'Drug Sale Arrests$_{t-1}$','Drug Possession Arrests$_{t-1}$','Domestic Arrests$_{t-1}$',
                  'Small Offenses Arrests$_{t-1}$','Missing Arrests$_{t-1}$',"Stock US items$_{t-1}$",
                  '$items/value * \frac{1}{dist_1}$', '$items/value * \frac{1}{dist_6}$', '$items/value * ln(landarea)$', '$items/value * 1[HIDTA]$',
                  'Stock US value$_{t-1}$','$items/value * \frac{1}{dist_1}$', '$items/value * \frac{1}{dist_6}$', '$items/value * ln(landarea)$', '$items/value * 1[HIDTA]$')

h.coef.order <- c(18:21, 1:17, 22)

texreg(harris.models,
       single.row = FALSE,
       custom.model.names = h.model.names,
       custom.coef.names = h.coef.names,
       reorder.coef = h.coef.order,
       caption = 'First Stage and Reduced Form Estimates of HPBM Results',
       label = "tab:reducedHarris",
       script.size=TRUE,
       caption.above = TRUE,
       file="tables/HPBM-HPBM-reduced-form.tex",
       ci.force=T,
       ci.test = NULL)

harris.county.models <- list(HPBM_first_stage_results[[3]], HPBM_county_red_items_results[[1]], HPBM_county_red_items_results[[2]],
                             HPBM_county_red_items_results[[3]], HPBM_county_red_items_results[[4]],
                             HPBM_first_stage_results[[4]], HPBM_county_red_value_results[[1]], HPBM_county_red_value_results[[2]],
                             HPBM_county_red_value_results[[3]], HPBM_county_red_value_results[[4]])

h.coef.names <- c("Percent Poverty","Unemployment Rate","log(Population)"," Murder Arrests$_{t-1}$",
                  "Rape Arrests$_{t-1}$","Robbery Arrests$_{t-1}$","Ag. Assault Arrests$_{t-1}$","Burglary Arrests$_{t-1}$",
                  "Assault Arrests$_{t-1}$","Weapons Arrests$_{t-1}$","Drug Sale Arrests$_{t-1}$","Drug Possession Arrests$_{t-1}$",
                  "Stock US items$_{t-1}$",
                  '$items/value * \frac{1}{dist_1}$', '$items/value * \frac{1}{dist_6}$', '$items/value * ln(landarea)$', '$items/value * 1[HIDTA]$',
                  'Stock US value$_{t-1}$','$items/value * \frac{1}{dist_1}$', '$items/value * \frac{1}{dist_6}$', '$items/value * ln(landarea)$', '$items/value * 1[HIDTA]$')

h.coef.order <- c(14:17,1:12,13,18)

texreg(harris.county.models,
       single.row = FALSE,
       custom.model.names = h.model.names,
       custom.coef.names = h.coef.names,
       reorder.coef = h.coef.order,
       caption = 'First Stage and Reduced Form Estimates of County Replication Results for HPBM',
       label = "tab:reducedCountyHarris",
       script.size=TRUE,
       caption.above = TRUE,
       file="tables/HPBM-county-reduced-form.tex",
       ci.force=T,
       ci.test = NULL)

harris.agency.models <- list(HPBM_first_stage_results[[1]], HPBM_agency_red_items_results[[1]], HPBM_agency_red_items_results[[2]],
                             HPBM_agency_red_items_results[[3]], HPBM_agency_red_items_results[[4]],
                             HPBM_first_stage_results[[2]], HPBM_agency_red_value_results[[1]], HPBM_agency_red_value_results[[2]],
                             HPBM_agency_red_value_results[[3]], HPBM_agency_red_value_results[[4]])

h.coef.names <- c("Percent Poverty","Unemployment Rate","log(Population)",
                  'Murder Arrests$_{t-1}$', 'Rape Arrests$_{t-1}$', "Robbery Arrests$_{t-1}$",
                  'Burglary Arrests$_{t-1}$', 'Assault Arrests$_{t-1}$', "Stock US items$_{t-1}$",
                  '$items/value * \frac{1}{dist_1}$', '$items/value * \frac{1}{dist_6}$', '$items/value * ln(landarea)$', '$items/value * 1[HIDTA]$',
                  'Stock US value$_{t-1}$','$items/value * \frac{1}{dist_1}$', '$items/value * \frac{1}{dist_6}$', '$items/value * ln(landarea)$', '$items/value * 1[HIDTA]$')

h.coef.order <- c(10:13,1:9,14)

texreg(harris.agency.models,
       single.row = FALSE,
       custom.model.names = h.model.names,
       custom.coef.names = h.coef.names,
       reorder.coef = h.coef.order,
       caption = 'First Stage and Reduced Form Estimates of Agency Replication Results for HPBM',
       label = "tab:reducedAgencyHarris",
       script.size=TRUE,
       caption.above = TRUE,
       file="tables/HPBM-agency-reduced-form.tex",
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 9 (Power)
## power calculation --------
load('data/results-BG.RData')
load('data/results-HPBM.RData')

# extract original estimates from BG ----
BG_orig_est <- unnest(BG_results, cols = c(`map(BG_list, tidy)`)) %>%
  filter(term == "`ltotal_cost(fit)`") %>%
  mutate(var = recode(term, "`ltotal_cost(fit)`" = "aid"),
         model = c('Total Crime Rate','Homicide', 'Robbery',
                   'Assault', 'Burglary','Larceny',"Vehicle Theft"),
         group = 'BG County') %>%
  select(estimate)

# extract standard error from our county estimates ---
BG_county_se <- unnest(BG_county_results, cols = c(`map(BG_county_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOaidBG_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidBG_lag)(fit)`" = "aid"),
         model = c('Total Crime Rate','Homicide', 'Robbery',
                   'Assault', 'Burglary','Larceny',"Vehicle Theft"),
         group = 'BG County Rep') %>%
  select(std.error)

# extract standard error from our agency estimates -----
BG_agency_se <- unnest(BG_agency_results, cols = c(agency)) %>%
  filter(term == "`log(1 + sumallLESOaidBG_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidBG_lag)(fit)`" = "aid"),
         model = c('Total Crime Rate','Homicide', 'Robbery',
                   'Assault', 'Burglary','Larceny',"Vehicle Theft"),
         group = 'BG Agency Rep') %>%
  select(std.error)

BG_orig_se <- unnest(BG_results, cols = c(`map(BG_list, tidy)`)) %>%
  filter(term == "`ltotal_cost(fit)`") %>%
  mutate(var = recode(term, "`ltotal_cost(fit)`" = "aid"),
         model = c('Total Crime Rate','Homicide', 'Robbery',
                   'Assault', 'Burglary','Larceny',"Vehicle Theft"),
         group = 'BG County') %>%
  select(std.error, p.value)

BG_est <- data.frame(beta = BG_orig_est$estimate,
                     se = BG_orig_se$std.error,
                     pvalue = BG_orig_se$p.value,
                     county_se = BG_county_se$std.error,
                     agency_se = BG_agency_se$std.error)

crime_types <- c("Total Crime", "Homicide", "Robbery", "Assault", "Burglary", "Larceny", "Vehicle Theft")

# get width of CI and estimate power
BG_est <- BG_est %>%
  mutate(se_2 = se*1.96,
         county_se_2 = county_se*1.96,
         agency_se_2 = agency_se*1.96,
         BG_power = pnorm(q = -se_2, mean = beta, sd = se),
         county_power = pnorm(q = -county_se_2, mean = beta, sd = county_se),
         agency_power = pnorm(q = -agency_se_2, mean = beta, sd = agency_se),
         Outcome = crime_types)

BG_table_data <- BG_est %>%
  select(Outcome,beta,pvalue,BG_power,county_power,agency_power) %>%
  mutate(beta=round(beta,2),
         orig_power = round(BG_power*100,1),
         county_power=round(county_power*100,1),
         agency_power=round(agency_power*100,1)) %>%
  select(-BG_power) %>%
  select(Outcome,beta,pvalue,orig_power,county_power, agency_power)

## Do HPBM data
# extract original main effect estimate ------
HPBM_items_est <- unnest(HPBM_items_results, cols = c(`map(harris_items_list, tidy)`)) %>%
  filter(term == "`log_sum_items_lag(fit)`") %>%
  mutate(var = recode(term, "`log_sum_items_lag(fit)`" = "sum_items"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Items') %>%
  select(estimate)

HPBM_value_est <- unnest(HPBM_value_results, cols = c(`map(harris_value_list, tidy)`)) %>%
  filter(term == "`log_sum_value_ttl_lag(fit)`") %>%
  mutate(var = recode(term, "`log_sum_value_ttl_lag(fit)`" = "sum_value"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Value') %>%
  select(estimate)

# extract our counte SE ------
HPBM_items_county_se <- unnest(HPBM_county_items_results, cols = c(`map(county_items_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "sum_items"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Rep Items') %>%
  select(std.error)

HPBM_value_county_se <- unnest(HPBM_county_value_results, cols = c(`map(county_value_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOaidHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "sum_value"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Rep Value') %>%
  select(std.error)

## extract our agency level SE ------
HPBM_items_agency_se <- unnest(HPBM_agency_items_results, cols = c(`map(agency_items_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "sum_items"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM Agency Rep Items') %>%
  select(std.error)

HPBM_value_agency_se <- unnest(HPBM_agency_value_results, cols = c(`map(agency_value_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOaidHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "sum_value"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM Agency Rep Value') %>%
  select(std.error)

HPBM_items_se <- unnest(HPBM_items_results, cols = c(`map(harris_items_list, tidy)`)) %>%
  filter(term == "`log_sum_items_lag(fit)`") %>%
  mutate(var = recode(term, "`log_sum_items_lag(fit)`" = "sum_items"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Items') %>%
  select(std.error, p.value)

HPBM_value_se <- unnest(HPBM_value_results, cols = c(`map(harris_value_list, tidy)`)) %>%
  filter(term == "`log_sum_value_ttl_lag(fit)`") %>%
  mutate(var = recode(term, "`log_sum_value_ttl_lag(fit)`" = "sum_value"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Value') %>%
  select(std.error, p.value)

HPBM_items_est <- data.frame(beta = HPBM_items_est$estimate,
                             se = HPBM_items_se$std.error,
                             pvalue = HPBM_items_se$p.value,
                             county_se = HPBM_items_county_se$std.error,
                             agency_se = HPBM_items_agency_se$std.error)

HPBM_value_est <- data.frame(beta = HPBM_value_est$estimate,
                             se = HPBM_value_se$std.error,
                             pvalue = HPBM_value_se$p.value,
                             county_se = HPBM_value_county_se$std.error,
                             agency_se = HPBM_value_agency_se$std.error)

crime_types <- c("Homicide", "Robbery", "Assault", "Vehicle Theft")

# get width of CI and estimate power
HPBM_items_est <- HPBM_items_est %>%
  mutate(se_2 = se*1.96,
         county_se_2 = county_se*1.96,
         agency_se_2 = agency_se*1.96,
         HPBM_power = pnorm(q = -se_2, mean = beta, sd = se),
         county_power = pnorm(q = -county_se_2, mean = beta, sd = county_se),
         agency_power = pnorm(q = -agency_se_2, mean = beta, sd = agency_se),
         Outcome = crime_types)

HPBM_value_est <- HPBM_value_est %>%
  mutate(se_2 = se*1.96,
         county_se_2 = county_se*1.96,
         agency_se_2 = agency_se*1.96,
         HPBM_power = pnorm(q = -se_2, mean = beta, sd = se),
         county_power = pnorm(q = -county_se_2, mean = beta, sd = county_se),
         agency_power = pnorm(q = -agency_se_2, mean = beta, sd = agency_se),
         Outcome = crime_types)

## fix homicide bc positive coefficient so CI should be other side
fix_orig <- 1 - pnorm(q = HPBM_value_est$se_2[1],mean = HPBM_value_est$beta[1],sd=HPBM_value_est$se[1])
fix_county <- 1 - pnorm(q = HPBM_value_est$county_se_2[1],mean = HPBM_value_est$beta[1],sd=HPBM_value_est$county_se[1])
fix_agency <- 1 - pnorm(q = HPBM_value_est$agency_se_2[1],mean = HPBM_value_est$beta[1],sd=HPBM_value_est$agency_se[1])
HPBM_value_est$HPBM_power[1] <- fix_orig
HPBM_value_est$county_power[1] <- fix_county
HPBM_value_est$agency_power[1] <- fix_agency

HPBM_table_data_items <- HPBM_items_est %>%
  select(Outcome,beta,pvalue,HPBM_power,county_power,agency_power) %>%
  mutate(beta=round(beta,2),
         orig_power = round(HPBM_power*100,1),
         county_power=round(county_power*100,1),
         agency_power=round(agency_power*100,1)) %>%
  select(-HPBM_power) %>%
  select(Outcome,beta,pvalue,orig_power,county_power, agency_power)

HPBM_table_data_value <- HPBM_value_est %>%
  select(Outcome,beta,pvalue,HPBM_power,county_power,agency_power) %>%
  mutate(beta=round(beta,2),
         orig_power = round(HPBM_power*100,1),
         county_power=round(county_power*100,1),
         agency_power=round(agency_power*100,1)) %>%
  select(-HPBM_power) %>%
  select(Outcome,beta,pvalue,orig_power,county_power, agency_power)

df <- bind_rows(BG_table_data, HPBM_table_data_items, HPBM_table_data_value)

print(xtable(df,caption="Power Analysis",label="tab:power_table"),
      math.style.exponents=FALSE,
      include.rownames=FALSE,
      type='latex',
      file='tables/power-analysis.tex')

rm(list=ls())


##Create Supplementary Table 10 (Equivalence Tests)
library(tidyverse)
library(xtable)

load('data/results-BG.RData')
load('data/results-HPBM.RData')


BG_orig_est <- unnest(BG_results, cols = c(`map(BG_list, tidy)`)) %>%
  filter(term == "`ltotal_cost(fit)`") %>%
  mutate(var = recode(term, "`ltotal_cost(fit)`" = "aid"),
         model = c('Total Crime Rate','Homicide', 'Robbery',
                   'Assault', 'Burglary','Larceny',"Vehicle Theft"),
         group = 'BG County') %>%
  select(estimate,std.error, p.value)

# extract standard error from our county estimates ---
BG_county_est <- unnest(BG_county_results, cols = c(`map(BG_county_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOaidBG_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidBG_lag)(fit)`" = "aid"),
         model = c('Total Crime Rate','Homicide', 'Robbery',
                   'Assault', 'Burglary','Larceny',"Vehicle Theft"),
         group = 'BG County Rep') %>%
  select(estimate,std.error)

# extract standard error from our agency estimates -----
BG_agency_est <- unnest(BG_agency_results, cols = c(agency)) %>%
  filter(term == "`log(1 + sumallLESOaidBG_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidBG_lag)(fit)`" = "aid"),
         model = c('Total Crime Rate','Homicide', 'Robbery',
                   'Assault', 'Burglary','Larceny',"Vehicle Theft"),
         group = 'BG Agency Rep') %>%
  select(estimate,std.error)


BG_est <- data.frame(beta = BG_orig_est$estimate,
                     se = BG_orig_est$std.error,
                     pvalue = BG_orig_est$p.value,
                     county_est = BG_county_est$estimate,
                     county_se = BG_county_est$std.error,
                     agency_est = BG_agency_est$estimate,
                     agency_se = BG_agency_est$std.error)

# get width of CI and estimate power
crime_types <- c("Total Crime", "Homicide", "Robbery", "Assault", "Burglary", "Larceny", "Vehicle Theft")

# 90% CI ------
alpha <- abs(qnorm(.05)) #90%

cty_low <- BG_est$county_est - alpha*BG_est$county_se
cty_high <- BG_est$county_est + alpha*BG_est$county_se
agcy_low <- BG_est$agency_est - alpha*BG_est$agency_se
agcy_high <- BG_est$agency_est + alpha*BG_est$agency_se


## HPBM data --------
# extract original main effect estimate HPBM ------
HPBM_items_est <- unnest(HPBM_items_results, cols = c(`map(harris_items_list, tidy)`)) %>%
  filter(term == "`log_sum_items_lag(fit)`") %>%
  mutate(var = recode(term, "`log_sum_items_lag(fit)`" = "sum_items"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Items') %>%
  select(estimate, std.error, p.value)

HPBM_value_est <- unnest(HPBM_value_results, cols = c(`map(harris_value_list, tidy)`)) %>%
  filter(term == "`log_sum_value_ttl_lag(fit)`") %>%
  mutate(var = recode(term, "`log_sum_value_ttl_lag(fit)`" = "sum_value"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Value') %>%
  select(estimate, std.error, p.value)

# extract our counte SE ------
HPBM_items_county_est <- unnest(HPBM_county_items_results, cols = c(`map(county_items_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "sum_items"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Rep Items') %>%
  select(estimate, std.error, p.value)

HPBM_value_county_est <- unnest(HPBM_county_value_results, cols = c(`map(county_value_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOaidHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "sum_value"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM County Rep Value') %>%
  select(estimate, std.error, p.value)

## extract our agency level SE ------
HPBM_items_agency_est <- unnest(HPBM_agency_items_results, cols = c(`map(agency_items_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "sum_items"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM Agency Rep Items') %>%
  select(estimate, std.error, p.value)

HPBM_value_agency_est <- unnest(HPBM_agency_value_results, cols = c(`map(agency_value_list, tidy)`)) %>%
  filter(term == "`log(1 + sumallLESOaidHARRIS_lag)(fit)`") %>%
  mutate(var = recode(term, "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "sum_value"),
         model = c('Homicide','Robbery','Assault',"Vehicle Theft"),
         group = 'HPBM Agency Rep Value') %>%
  select(estimate, std.error, p.value)

pvalue <- HPBM_items_est$p.value
HPBM_items <- data.frame(beta = HPBM_items_est$estimate,
                         se = HPBM_items_est$std.error,
                         county_est = HPBM_items_county_est$estimate,
                         cty_low = HPBM_items_county_est$estimate - alpha*HPBM_items_county_est$std.error,
                         cty_high = HPBM_items_county_est$estimate + alpha*HPBM_items_county_est$std.error,
                         agency_est = HPBM_items_agency_est$estimate,
                         agcy_low = HPBM_items_agency_est$estimate - alpha*HPBM_items_agency_est$std.error,
                         agcy_high = HPBM_items_agency_est$estimate + alpha*HPBM_items_agency_est$std.error) %>%
  round(.,2) %>%
  bind_cols(.,data.frame(pvalue)) %>%
  select(beta,se,pvalue,everything())

pvalue <- HPBM_value_est$p.value
HPBM_value <- data.frame(beta = HPBM_value_est$estimate,
                         se = HPBM_value_est$std.error,
                         county_est = HPBM_value_county_est$estimate,
                         cty_low = HPBM_value_county_est$estimate - alpha*HPBM_value_county_est$std.error,
                         cty_high = HPBM_value_county_est$estimate + alpha*HPBM_value_county_est$std.error,
                         agency_est = HPBM_value_agency_est$estimate,
                         agcy_low = HPBM_value_agency_est$estimate - alpha*HPBM_value_agency_est$std.error,
                         agcy_high = HPBM_value_agency_est$estimate + alpha*HPBM_value_agency_est$std.error) %>%
  round(.,2) %>%
  bind_cols(.,data.frame(pvalue)) %>%
  select(beta,se,pvalue,everything())

## table
pvalue <- BG_est$pvalue
BG_table_data <- BG_est %>%
  select(beta,se,county_est,agency_est) %>%
  bind_cols(data.frame(cty_low,cty_high,agcy_low,agcy_high)) %>%
  round(.,2) %>%
  bind_cols(data.frame(pvalue)) %>%
  select(beta,se,pvalue,county_est,cty_low,cty_high,agency_est,agcy_low,agcy_high)


full_table <- bind_rows(BG_table_data, HPBM_items, HPBM_value)

full_table <- full_table %>% select(-se,-cty_high, -agcy_high)
outcome <- data.frame(Outcome = c(crime_types,rep(c("Homicide",'Robbery','Assault','Vehicle Theft'),2)))
full_table <- bind_cols(outcome,full_table)

print(xtable(full_table, display = c("s","s","g","g","s","s","s","s"),
             caption="Equivalence Test",label="tab:equiv_test"),
      math.style.exponents=FALSE,
      include.rownames=FALSE,
      type='latex',
      file='tables/CI_90percent.tex')

rm(list=ls())


## Create Supplementary Table 11 (BG Controlled Items)
load('data/results-BG-for-tables-controlled.RData')

table2.A <- list(BG_OLS_results[[3]], BG_OLS_results[[2]],BG_OLS_results[[1]],
                 first_stage_results[[1]], first_stage_results[[3]], first_stage_results[[2]],
                 BG_results[[1]], BG_county_results[[1]], BG_agency_results[[1]])

table2.B <- list(BG_results[[2]], BG_county_results[[2]], BG_agency_results[[2]],
                 BG_results[[3]], BG_county_results[[3]], BG_agency_results[[3]],
                 BG_results[[4]], BG_county_results[[4]], BG_agency_results[[4]])

table2.C <- list(BG_results[[5]], BG_county_results[[5]], BG_agency_results[[5]],
                 BG_results[[6]], BG_county_results[[6]], BG_agency_results[[6]],
                 BG_results[[7]], BG_county_results[[7]], BG_agency_results[[7]])

coef.map <- list(milex_iv = "Military exp. IV",
                 ltotal_cost = "Lagged Total Aid",
                 "log(1 + sumallLESOaidBG_lag)" = "Lagged Total Aid",
                 "`log(1 + sumallLESOaidBG_lag)(fit)`" = "Lagged Total Aid",
                 "`ltotal_cost(fit)`" = "Lagged Total Aid")

texreg(table2.A,
       file = "tables/BGtable2-partA-new-controlled.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.B,
       file = "tables/BGtable2-partB-new-controlled.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.C,
       file = "tables/BGtable2-partC-new-controlled.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 12 (HPBM Controlled Items)
load('data/results-HPBM-for-tables-controlled.RData')

partA.1 <- list(HPBM_items_results[[1]], HPBM_county_items_results[[1]], HPBM_agency_items_results[[1]],
                HPBM_items_results[[2]], HPBM_county_items_results[[2]], HPBM_agency_items_results[[2]])

partA.2 <- list(HPBM_value_results[[1]], HPBM_county_value_results[[1]], HPBM_agency_value_results[[1]],
                HPBM_value_results[[2]], HPBM_county_value_results[[2]], HPBM_agency_value_results[[2]])

partB.1 <- list(HPBM_items_results[[3]], HPBM_county_items_results[[3]], HPBM_agency_items_results[[3]],
                HPBM_items_results[[4]], HPBM_county_items_results[[4]], HPBM_agency_items_results[[4]])

partB.2 <- list(HPBM_value_results[[3]], HPBM_county_value_results[[3]], HPBM_agency_value_results[[3]],
                HPBM_value_results[[4]], HPBM_county_value_results[[4]], HPBM_agency_value_results[[4]])

coef.map1 <- list("`log_sum_items_lag(fit)`" = "log items$_{t-1}$",
                  "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "log items$_{t-1}$")
coef.map2 <- list("`log_sum_value_ttl_lag(fit)`" = "log value$_{t-1}$",
                  "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "log value$_{t-1}$")

texreg(partA.1,
       file = "tables/HPBMtable8-partA1-controlled.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partA.2,
       file = "tables/HPBMtable8-partA2-controlled.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.1,
       file = "tables/HPBMtable8-partB1-controlled.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.2,
       file = "tables/HPBMtable8-partB2-controlled.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 13 (Main BG Results 2010-2012 Subset)
load('data/results-BG-for-tables-year-subset.RData')

table2.A <- list(BG_OLS_results[[3]], BG_OLS_results[[2]],BG_OLS_results[[1]],
                 first_stage_results[[1]], first_stage_results[[3]], first_stage_results[[2]],
                 BG_results[[1]], BG_county_results[[1]], BG_agency_results[[1]])

table2.B <- list(BG_results[[2]], BG_county_results[[2]], BG_agency_results[[2]],
                 BG_results[[3]], BG_county_results[[3]], BG_agency_results[[3]],
                 BG_results[[4]], BG_county_results[[4]], BG_agency_results[[4]])

table2.C <- list(BG_results[[5]], BG_county_results[[5]], BG_agency_results[[5]],
                 BG_results[[6]], BG_county_results[[6]], BG_agency_results[[6]],
                 BG_results[[7]], BG_county_results[[7]], BG_agency_results[[7]])

coef.map <- list(milex_iv = "Military exp. IV",
                 ltotal_cost = "Lagged Total Aid",
                 "log(1 + sumallLESOaidBG_lag)" = "Lagged Total Aid",
                 "`log(1 + sumallLESOaidBG_lag)(fit)`" = "Lagged Total Aid",
                 "`ltotal_cost(fit)`" = "Lagged Total Aid")

texreg(table2.A,
       file = "tables/BGtable2-partA-new-subset.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.B,
       file = "tables/BGtable2-partB-new-subset.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.C,
       file = "tables/BGtable2-partC-new-subset.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 14 (Main HPBM Results 2010-2012 Subset)
load('data/results-HPBM-for-tables-year-subset.RData')

partA.1 <- list(HPBM_items_results[[1]], HPBM_county_items_results[[1]], HPBM_agency_items_results[[1]],
                HPBM_items_results[[2]], HPBM_county_items_results[[2]], HPBM_agency_items_results[[2]])

partA.2 <- list(HPBM_value_results[[1]], HPBM_county_value_results[[1]], HPBM_agency_value_results[[1]],
                HPBM_value_results[[2]], HPBM_county_value_results[[2]], HPBM_agency_value_results[[2]])

partB.1 <- list(HPBM_items_results[[3]], HPBM_county_items_results[[3]], HPBM_agency_items_results[[3]],
                HPBM_items_results[[4]], HPBM_county_items_results[[4]], HPBM_agency_items_results[[4]])

partB.2 <- list(HPBM_value_results[[3]], HPBM_county_value_results[[3]], HPBM_agency_value_results[[3]],
                HPBM_value_results[[4]], HPBM_county_value_results[[4]], HPBM_agency_value_results[[4]])

coef.map1 <- list("`log_sum_items_lag(fit)`" = "log items$_{t-1}$",
                  "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "log items$_{t-1}$")
coef.map2 <- list("`log_sum_value_ttl_lag(fit)`" = "log value$_{t-1}$",
                  "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "log value$_{t-1}$")

texreg(partA.1,
       file = "tables/HPBMtable8-partA1-subset.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partA.2,
       file = "tables/HPBMtable8-partA2-subset.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.1,
       file = "tables/HPBMtable8-partB1-subset.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.2,
       file = "tables/HPBMtable8-partB2-subset.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 15 (Main BG Results No Outliers Elminated)
load('data/results-BG-for-tables-with-outliers.RData')

table2.A <- list(BG_OLS_results[[3]], BG_OLS_results[[2]],BG_OLS_results[[1]],
                 first_stage_results[[1]], first_stage_results[[3]], first_stage_results[[2]],
                 BG_results[[1]], BG_county_results[[1]], BG_agency_results[[1]])

table2.B <- list(BG_results[[2]], BG_county_results[[2]], BG_agency_results[[2]],
                 BG_results[[3]], BG_county_results[[3]], BG_agency_results[[3]],
                 BG_results[[4]], BG_county_results[[4]], BG_agency_results[[4]])

table2.C <- list(BG_results[[5]], BG_county_results[[5]], BG_agency_results[[5]],
                 BG_results[[6]], BG_county_results[[6]], BG_agency_results[[6]],
                 BG_results[[7]], BG_county_results[[7]], BG_agency_results[[7]])

coef.map <- list(milex_iv = "Military exp. IV",
                 ltotal_cost = "Lagged Total Aid",
                 "log(1 + sumallLESOaidBG_lag)" = "Lagged Total Aid",
                 "`log(1 + sumallLESOaidBG_lag)(fit)`" = "Lagged Total Aid",
                 "`ltotal_cost(fit)`" = "Lagged Total Aid")

texreg(table2.A,
       file = "tables/BGtable2-partA-new-outliers.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.B,
       file = "tables/BGtable2-partB-new-outliers.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.C,
       file = "tables/BGtable2-partC-new-outliers.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 16 (Main HPBM Results No Outliers Eliminated)
load('data/results-HPBM-for-tables-with-outliers.RData')

partA.1 <- list(HPBM_items_results[[1]], HPBM_county_items_results[[1]], HPBM_agency_items_results[[1]],
                HPBM_items_results[[2]], HPBM_county_items_results[[2]], HPBM_agency_items_results[[2]])

partA.2 <- list(HPBM_value_results[[1]], HPBM_county_value_results[[1]], HPBM_agency_value_results[[1]],
                HPBM_value_results[[2]], HPBM_county_value_results[[2]], HPBM_agency_value_results[[2]])

partB.1 <- list(HPBM_items_results[[3]], HPBM_county_items_results[[3]], HPBM_agency_items_results[[3]],
                HPBM_items_results[[4]], HPBM_county_items_results[[4]], HPBM_agency_items_results[[4]])

partB.2 <- list(HPBM_value_results[[3]], HPBM_county_value_results[[3]], HPBM_agency_value_results[[3]],
                HPBM_value_results[[4]], HPBM_county_value_results[[4]], HPBM_agency_value_results[[4]])

coef.map1 <- list("`log_sum_items_lag(fit)`" = "log items$_{t-1}$",
                  "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "log items$_{t-1}$")
coef.map2 <- list("`log_sum_value_ttl_lag(fit)`" = "log value$_{t-1}$",
                  "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "log value$_{t-1}$")

texreg(partA.1,
       file = "tables/HPBMtable8-partA1-outliers.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partA.2,
       file = "tables/HPBMtable8-partA2-outliers.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.1,
       file = "tables/HPBMtable8-partB1-outliers.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.2,
       file = "tables/HPBMtable8-partB2-outliers.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 17 (Main BG Results Dropping Problematic Crime Observations)
load('data/results-BG-for-tables-crime-obs-dropped.RData')

table2.A <- list(BG_OLS_results[[3]], BG_OLS_results[[2]],BG_OLS_results[[1]],
                 first_stage_results[[1]], first_stage_results[[3]], first_stage_results[[2]],
                 BG_results[[1]], BG_county_results[[1]], BG_agency_results[[1]])

table2.B <- list(BG_results[[2]], BG_county_results[[2]], BG_agency_results[[2]],
                 BG_results[[3]], BG_county_results[[3]], BG_agency_results[[3]],
                 BG_results[[4]], BG_county_results[[4]], BG_agency_results[[4]])

table2.C <- list(BG_results[[5]], BG_county_results[[5]], BG_agency_results[[5]],
                 BG_results[[6]], BG_county_results[[6]], BG_agency_results[[6]],
                 BG_results[[7]], BG_county_results[[7]], BG_agency_results[[7]])

coef.map <- list(milex_iv = "Military exp. IV",
                 ltotal_cost = "Lagged Total Aid",
                 "log(1 + sumallLESOaidBG_lag)" = "Lagged Total Aid",
                 "`log(1 + sumallLESOaidBG_lag)(fit)`" = "Lagged Total Aid",
                 "`ltotal_cost(fit)`" = "Lagged Total Aid")

texreg(table2.A,
       file = "tables/BGtable2-partA-new-crime-obs-dropped.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.B,
       file = "tables/BGtable2-partB-new-crime-obs-dropped.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(table2.C,
       file = "tables/BGtable2-partC-new-crime-obs-dropped.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = rep(c("BG","County","Agency"),3),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 18 (Main HPBM Results Dropping Problematic Crime Observations)
load('data/results-HPBM-for-tables-crime-obs-dropped.RData')

partA.1 <- list(HPBM_items_results[[1]], HPBM_county_items_results[[1]], HPBM_agency_items_results[[1]],
                HPBM_items_results[[2]], HPBM_county_items_results[[2]], HPBM_agency_items_results[[2]])

partA.2 <- list(HPBM_value_results[[1]], HPBM_county_value_results[[1]], HPBM_agency_value_results[[1]],
                HPBM_value_results[[2]], HPBM_county_value_results[[2]], HPBM_agency_value_results[[2]])

partB.1 <- list(HPBM_items_results[[3]], HPBM_county_items_results[[3]], HPBM_agency_items_results[[3]],
                HPBM_items_results[[4]], HPBM_county_items_results[[4]], HPBM_agency_items_results[[4]])

partB.2 <- list(HPBM_value_results[[3]], HPBM_county_value_results[[3]], HPBM_agency_value_results[[3]],
                HPBM_value_results[[4]], HPBM_county_value_results[[4]], HPBM_agency_value_results[[4]])

coef.map1 <- list("`log_sum_items_lag(fit)`" = "log items$_{t-1}$",
                  "`log(1 + sumallLESOquantityHARRIS_lag)(fit)`" = "log items$_{t-1}$")
coef.map2 <- list("`log_sum_value_ttl_lag(fit)`" = "log value$_{t-1}$",
                  "`log(1 + sumallLESOaidHARRIS_lag)(fit)`" = "log value$_{t-1}$")

texreg(partA.1,
       file = "tables/HPBMtable8-partA1-crime-obs-dropped.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partA.2,
       file = "tables/HPBMtable8-partA2-crime-obs-dropped.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.1,
       file = "tables/HPBMtable8-partB1-crime-obs-dropped.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.2,
       file = "tables/HPBMtable8-partB2-crime-obs-dropped.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("HPBM", "County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

#Clear working environment
rm(list=ls())


## Create Supplementary Table 19 (Alternative Approach to Calculating Quantity and Value for HPBM)
# run accumulated sum for HPBM --------
# cleaned data
library(plm) #for panel data
library(lfe) #for felm function
library(rlang)
library(tidyverse)
library(broom)
library(foreign)
library(purrr)
library(texreg)
library(foreign)
library(readxl)
library(openintro)
library(lubridate)
source("lib/model-functions-HPBM.R") #these have all the different models we need to run
source("lib/fix-names.R")

load('data/city-aggregate.RData')
load('data/county-aggregate.RData')
city_outliers_eliminated <- TRUE
if(city_outliers_eliminated){
  agency <- agency %>% filter(allcrime_rate < 1000000)
}

#----------------------------------------------------------------#
# Military Aid Data 
#----------------------------------------------------------------#
path <- "raw-Data/DISP_AllStatesAndTerritories_03312018.xlsx"
LESOlist <- path %>% 
  excel_sheets() %>% 
  set_names() %>% 
  map(read_excel, path = path)

LESO <- do.call(rbind,LESOlist) %>% 
  rename(agency_name=`Station Name (LEA)`) %>%
  mutate(state_abbrv = State,
         State=abbr2state(State),
         `Ship Date` =  ymd_hms(`Ship Date`),
         year = year(`Ship Date`)) %>%
  filter(!is.na(State))

LESO$total_value <- LESO$Quantity * LESO$`Acquisition Value`

# aggregate to the year by agency
leso_year_sum <- LESO %>%
  group_by(State, agency_name, year) %>%
  summarize(total_year_value = sum(total_value),
            total_year_quantity = sum(Quantity)) %>%
  ungroup() %>%
  group_by(State,agency_name)


# need to merge in city and county info so I can merge with 
# already cleaned military dataset


# Read in fips data and match by agency name 
fips <- read_csv("raw-data/clean-countiesandfips.csv") %>%
  select(agency,agency_name,State,city,County)

#militarization
city <- left_join(leso_year_sum,fips,
                  by = c("State"="State",
                         "agency_name"="agency_name")) %>%
  fix_county_names(.) %>%
  select(-agency)

# EC "NAVAJO POLICE DEPT, TUBA CTY DIST" is in military data but not counties
# they had 3 trucks
city$city[city$agency_name=="NAVAJO POLICE DEPT, TUBA CTY DIST"] <- "Tuba City"
city$County[city$agency_name=="NAVAJO POLICE DEPT, TUBA CTY DIST"] <- "Coconino County"

city$city[city$agency_name=="WI DEPT OF JUSTICE, DCI (LEA)"] <- "Madison"
city$County[city$agency_name=="WI DEPT OF JUSTICE, DCI (LEA)"] <- "Dane County"

# create all year combinations
# start from earliest year bc need stock
city_year_combo <- city %>%
  tidyr::expand(State, city, 1990:2015) %>%
  rename(year = `1990:2015`)

city_year <- left_join(city_year_combo, city) %>%
  select(-County)

city_year$total_year_value[is.na(city_year$total_year_value)] <- 0
city_year$total_year_quantity[is.na(city_year$total_year_quantity)] <- 0
# need aggregated to city 
city_year_sum <- city_year %>%
  group_by(State, city, year) %>%
  summarize(total_year_value = sum(total_year_value),
            total_year_quantity = sum(total_year_quantity)) %>%
  mutate(cum_value = cumsum(total_year_value),
         cum_quantity = cumsum(total_year_quantity),
         cum_value_lag = lag(cum_value),
         cum_quantity_lag = lag(cum_quantity)) %>%
  filter(year >= 2010 & year <= 2015)


agency <- agency %>%
  select(City, State, year, PCPI, unemployment, Population, Murder_ArrestRate_lag,
         Rape_ArrestRate_lag, Robbery_ArrestRate_lag, Burglary_ArrestRate_lag, Assault_ArrestRate_lag,
         g, sumallLESOquantityHARRIS_lag, SumUSitems_lag, zHUdI1_lag, zHUdIland_lag,
         zHUdI6_lag, zHUdIhidta_lag, sumallLESOaidHARRIS_lag, SumUSvalue_lag,
         zHUd1_lag, zHUd6_lag, zHUdland_lag, zHUdhitda_lag,
         Murder_rate, Robbery_rate, Assault_rate, MV_theft_rate)

agency <- left_join(agency,city_year_sum,
                    by=c("City" = "city",
                         "State" = "State",
                         "year" = "year"))
agency$cum_value_lag[is.na(agency$cum_value_lag)] <- 0
agency$cum_quantity_lag[is.na(agency$cum_quantity_lag)] <- 0

# need aggregated to county 
# 1. aggregate to county
county_year <- city %>%
  group_by(State,County,year) %>%
  summarize(total_year_value = sum(total_year_value),
            total_year_quantity = sum(total_year_quantity))

#expand for all years
county_year_combo <- county_year %>%
  tidyr::expand(State, County, 1990:2015) %>%
  rename(year = `1990:2015`)

county_year <- left_join(county_year_combo, county_year)

county_year$total_year_value[is.na(county_year$total_year_value)] <- 0
county_year$total_year_quantity[is.na(county_year$total_year_quantity)] <- 0
# calculate cumulative stock county 
county_year_sum <- county_year %>%
  group_by(State, County) %>%
  mutate(cum_value = cumsum(total_year_value),
         cum_quantity = cumsum(total_year_quantity),
         cum_value_lag = lag(cum_value),
         cum_quantity_lag = lag(cum_quantity)) %>%
  filter(year >= 2010 & year <= 2015)

county <- county %>%
  select(County, State, year, PCPI,unemployment,Population, MURDER_ArrestRate_lag,
         RAPE_ArrestRate_lag,ROBBERY_ArrestRate_lag,AGASSLT_ArrestRate_lag,
         BURGLRY_ArrestRate_lag,OTHASLT_ArrestRate_lag,WEAPONS_ArrestRate_lag,
         DRGSALE_ArrestRate_lag,DRGPOSS_ArrestRate_lag, g,
         SumUSitems_lag, zHUdI1_lag, zHUdIland_lag, zHUdI6_lag, zHUdIhidta_lag,
         SumUSvalue_lag, zHUd1_lag, zHUd6_lag, zHUdland_lag, zHUdhitda_lag,
         MURDER_countysumrate, ROBBERY_countysumrate,
         AGASSLT_countysumrate, MVTHEFT_countysumrate,
         sumallLESOaidHARRIS_lag,sumallLESOquantityHARRIS_lag)

county <- left_join(county, county_year_sum,
                    by=c("County" = "County",
                         "State" = "State",
                         "year" = "year"))
county$cum_quantity_lag[is.na(county$cum_quantity_lag)] <- 0
county$cum_value_lag[is.na(county$cum_value_lag)] <- 0
## functions used in running models HPBM
## 2SLS 
run.IV.agency.Harris.cum <- function(outcome,data,type){
  outcome <- enexpr(outcome)
  if(type=="items"){
    form <-  as.formula(paste(quo_text(outcome),
                              '~ PCPI+ unemployment+Population+Murder_ArrestRate_lag +
                              Rape_ArrestRate_lag+Robbery_ArrestRate_lag+Burglary_ArrestRate_lag+
                              Assault_ArrestRate_lag', '| g+year',
                              '|(log(1+cum_quantity_lag)~SumUSitems_lag+zHUdI1_lag+ 
                              zHUdIland_lag+zHUdI6_lag+zHUdIhidta_lag)|g'))
  }else if(type=="value"){
    form <-  as.formula(paste(quo_text(outcome),
                              '~ PCPI+ unemployment+Population+Murder_ArrestRate_lag +
                              Rape_ArrestRate_lag+Robbery_ArrestRate_lag+Burglary_ArrestRate_lag+
                              Assault_ArrestRate_lag',
                              '| g+year',
                              '|(log(1+cum_value_lag)~SumUSvalue_lag+zHUd1_lag+zHUd6_lag+
                              zHUdland_lag+zHUdhitda_lag)|g'))
  }
  
  data %>% felm(form,data=.)
  }


## County level 
run.IV.county.Harris.cum <- function(outcome,data,type){
  outcome <- enexpr(outcome)
  county_vars_all_regressions <- paste(c("PCPI","unemployment","Population", "MURDER_ArrestRate_lag",
                                         "RAPE_ArrestRate_lag","ROBBERY_ArrestRate_lag","AGASSLT_ArrestRate_lag",
                                         "BURGLRY_ArrestRate_lag","OTHASLT_ArrestRate_lag","WEAPONS_ArrestRate_lag",
                                         "DRGSALE_ArrestRate_lag","DRGPOSS_ArrestRate_lag"),collapse='+')
  
  fe_variables_items <- paste('| g+year',
                              '|(log(1+cum_quantity_lag)~SumUSitems_lag+zHUdI1_lag+ 
                              zHUdIland_lag+zHUdI6_lag+zHUdIhidta_lag)|g')
  fe_variables_value <- paste('| g+year',
                              '|(log(1+cum_value_lag)~SumUSvalue_lag+zHUd1_lag+zHUd6_lag+
                              zHUdland_lag+zHUdhitda_lag)|g')
  if(type=="items"){
    form <-  as.formula(paste(quo_text(outcome),'~',county_vars_all_regressions,fe_variables_items))
  }else if(type=="value"){
    form <-  as.formula(paste(quo_text(outcome),'~',county_vars_all_regressions,fe_variables_value))
  }
  
  data %>% felm(form,data=.)
}

# create panel data
agency <- pdata.frame(agency,index=c("g","year"))
county <- pdata.frame(county,index=c("g","year"))
## Harris Agency level using cumulative 
homicide.items.agency.Harris <- run.IV.agency.Harris.cum(Murder_rate, agency, 'items')
homicide.value.agency.Harris <- run.IV.agency.Harris.cum(Murder_rate, agency, 'value')
robbery.items.agency.Harris <- run.IV.agency.Harris.cum(Robbery_rate, agency, 'items')
robbery.value.agency.Harris <- run.IV.agency.Harris.cum(Robbery_rate, agency, 'value')
assault.items.agency.Harris <- run.IV.agency.Harris.cum(Assault_rate, agency, 'items')
assault.value.agency.Harris <- run.IV.agency.Harris.cum(Assault_rate, agency, 'value')
vehicle.items.agency.Harris <- run.IV.agency.Harris.cum(MV_theft_rate, agency, 'items')
vehicle.value.agency.Harris <- run.IV.agency.Harris.cum(MV_theft_rate, agency, 'value')

## Harris County level using cumulative 
homicide.county.items.harris <- run.IV.county.Harris.cum(MURDER_countysumrate, county, 'items')
homicide.county.value.harris <- run.IV.county.Harris.cum(MURDER_countysumrate, county, 'value')
robbery.county.items.harris <- run.IV.county.Harris.cum(ROBBERY_countysumrate, county, 'items')
robbery.county.value.harris <- run.IV.county.Harris.cum(ROBBERY_countysumrate, county, 'value')
assault.county.items.harris <- run.IV.county.Harris.cum(AGASSLT_countysumrate, county, 'items')
assault.county.value.harris <- run.IV.county.Harris.cum(AGASSLT_countysumrate, county, 'value')
vehicle.county.items.harris <- run.IV.county.Harris.cum(MVTHEFT_countysumrate, county, 'items')
vehicle.county.value.harris <- run.IV.county.Harris.cum(MVTHEFT_countysumrate, county, 'value')


#---------------------------------------------------------#
## First stage 
#---------------------------------------------------------#
harris_agency_items_first <- summary(homicide.items.agency.Harris$stage1)$coefficients
harris_agency_value_first <- summary(homicide.value.agency.Harris$stage1)$coefficients
harris_county_items_first <- summary(homicide.county.items.harris$stage1)$coefficients
harris_county_value_first <- summary(homicide.county.value.harris$stage1)$coefficients

first_list <- list(harris_agency_items_first,
                   harris_agency_value_first,
                   harris_county_items_first,
                   harris_county_value_first)
HPBM_first_stage_results <- tibble(first_stage = map(first_list, tidy))

# Save all results 
agency_items_list <- list(homicide.items.agency.Harris,
                          robbery.items.agency.Harris,
                          assault.items.agency.Harris,
                          vehicle.items.agency.Harris)
agency_value_list <- list(homicide.value.agency.Harris,
                          robbery.value.agency.Harris,
                          assault.value.agency.Harris,
                          vehicle.value.agency.Harris)

county_items_list <- list(homicide.county.items.harris,
                          robbery.county.items.harris,
                          assault.county.items.harris,
                          vehicle.county.items.harris)
county_value_list <- list(homicide.county.value.harris,
                          robbery.county.value.harris,
                          assault.county.value.harris,
                          vehicle.county.value.harris)


HPBM_agency_items_results <- tibble(map(agency_items_list, tidy))
HPBM_agency_value_results <- tibble(map(agency_value_list, tidy))
HPBM_county_items_results <- tibble(map(county_items_list, tidy))
HPBM_county_value_results <- tibble(map(county_value_list, tidy))


# Save all results for texreg
HPBM_agency_items_results <- map(agency_items_list, texreg::extract)
HPBM_agency_value_results <- map(agency_value_list, texreg::extract)
HPBM_county_items_results <- map(county_items_list, texreg::extract)
HPBM_county_value_results <- map(county_value_list, texreg::extract)


harris_agency_items_first <- homicide.items.agency.Harris$stage1
harris_agency_value_first <- homicide.value.agency.Harris$stage1
harris_county_items_first <- homicide.county.items.harris$stage1
harris_county_value_first <- homicide.county.value.harris$stage1

first_list <- list(harris_agency_items_first,
                   harris_agency_value_first,
                   harris_county_items_first,
                   harris_county_value_first)
HPBM_first_stage_results <- map(first_list, texreg::extract)

partA.1 <- list(HPBM_county_items_results[[1]], HPBM_agency_items_results[[1]],
                HPBM_county_items_results[[2]], HPBM_agency_items_results[[2]])

partA.2 <- list(HPBM_county_value_results[[1]], HPBM_agency_value_results[[1]],
                HPBM_county_value_results[[2]], HPBM_agency_value_results[[2]])

partB.1 <- list(HPBM_county_items_results[[3]], HPBM_agency_items_results[[3]],
                HPBM_county_items_results[[4]], HPBM_agency_items_results[[4]])

partB.2 <- list(HPBM_county_value_results[[3]], HPBM_agency_value_results[[3]],
                HPBM_county_value_results[[4]], HPBM_agency_value_results[[4]])

coef.map1 <- list("`log(1 + cum_quantity_lag)(fit)`" = "log stock items$_{t-1}$")
coef.map2 <- list("`log(1 + cum_value_lag)(fit)`" = "log value$_{t-1}$")

texreg(partA.1,
       file = "tables/HPBMtable8-partA1-cum.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partA.2,
       file = "tables/HPBMtable8-partA2-cum.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.1,
       file = "tables/HPBMtable8-partB1-cum.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map1,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

texreg(partB.2,
       file = "tables/HPBMtable8-partB2-cum.tex",
       use.packages = FALSE,
       custom.coef.map = coef.map2,
       digits = 2,
       caption.above = TRUE,
       caption = "The Effect of Receiving Tactical Items on Substantiated Crime Rates",
       custom.model.names = rep(c("County","Agency"),2),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       ci.force=T,
       ci.test = NULL)

rm(list = ls())


## Create Supplementary Table 20 (Summary Statistics)
#Load in the data
load('data/BG_milit.RData')
load('data/harris.RData')
load('data/city-aggregate.RData')
load('data/county-aggregate.RData')
harris <- df

#Summary stats function
my.summary = function(x, na.rm=TRUE){
  result = c(Mean=mean(x, na.rm=na.rm),
             SD=sd(x, na.rm=na.rm),
             Median=median(x, na.rm=na.rm))
}

#Getting the summary stats for BG
BG_milit <- mutate(BG_milit, shareage2534 = shareage2529+shareage3034)
bg_ind <- c("crime", "murder", "robbery", "assault", "burglary", 
            "larceny", "mvtheft", "tot_q", "total_raw",           
            "PovertyPercentAllAges", "MedianHouseholdIncome", "Unempl_Rate", "pop", "sharemale",            
            "shareblack", "shareage1519", "shareage2024", "shareage2534", "milex_iv")
bg_summary <- as.data.frame(sapply(BG_milit[, bg_ind], my.summary))

#Getting the summary stats for Harris
harris <- mutate(harris, sum_items = exp(log_sum_items))
harris <- mutate(harris, sum_items = if_else(sum_items==1, 0, sum_items))
harris <- mutate(harris, sum_value = exp(log_sum_value_ttl))
harris <- mutate(harris, sum_value = if_else(sum_value==1, 0, sum_value))
harris_ind <- c("Amurders", "Arobbery", "Aassault", "Aburglar", "Avehicle", "PCPI", 
                "UnemploymentRate", "allpop", "blackshare", "sum_value", "sum_items", "zHUdI1",           
                "zHUdI6", "zHUdIland", "zHUdIhidta", "zHUd1", "zHUd6",            
                "zHUdland", "zHUdhitda")
harris_summary <- as.data.frame(sapply(harris[, harris_ind], my.summary))

#Getting the summary stats for Emory Jurisdiction (Agency)
agency <- mutate(agency, shareblack = percblack/100)
agency_ind <- c("allcrime_rate", "Murder_rate", "Robbery_rate", "Assault_rate", "Burglary_rate",           
                "Larceny_rate", "MV_theft_rate", "sumallLESOaidBG", "sumallLESOaidHARRIS",     
                "sumallLESOquantityBG", "sumallLESOquantityHARRIS", "percentinpoverty", "medianhouseholdincome", "unemployment",             
                "Total_Population", "sharemales", "shareblack", "share1519", "share2024",             
                "share2534", "PCPI", "milex_iv", "zHUdI1",           
                "zHUdI6", "zHUdIland", "zHUdIhidta", "zHUd1", "zHUd6",            
                "zHUdland", "zHUdhitda")
agency_summary <- as.data.frame(sapply(agency[, agency_ind], my.summary))

#Getting the summary stats for Emory County
county <- mutate(county, shareblack = percblack/100)
county_ind <- c("allcrime_countysumrate", "MURDER_countysumrate", "ROBBERY_countysumrate", "AGASSLT_countysumrate", "BURGLRY_countysumrate",        
                "LARCENY_countysumrate", "MVTHEFT_countysumrate", "sumallLESOaidBG", "sumallLESOaidHARRIS",      
                "sumallLESOquantityBG", "sumallLESOquantityHARRIS", "percentinpoverty", "medianhouseholdincome", "unemployment",                 
                "Population", "sharemales", "shareblack", "share1519", "share2024",                 
                "share2534", "PCPI", "milex_iv", "zHUdI1",           
                "zHUdI6", "zHUdIland", "zHUdIhidta", "zHUd1", "zHUd6",            
                "zHUdland", "zHUdhitda")
county_summary <- as.data.frame(sapply(county[, county_ind], my.summary))

#Making the table
agency_means <- c(agency_summary[1,])
agency_sds <- c(agency_summary[2,])
agency_medians <- c(agency_summary[3,])
county_means <- c(county_summary[1,])
county_sds <- c(county_summary[2,])
county_medians <- c(county_summary[3,])
bg_means <- c(bg_summary[1,1:7], bg_summary[1,9], NA, bg_summary[1,8], NA, bg_summary[1, 10:18], NA, bg_summary[1, 19], rep(NA, 8))
bg_sds <- c(bg_summary[2,1:7], bg_summary[2,9], NA, bg_summary[2,8], NA, bg_summary[2, 10:18], NA, bg_summary[2,19], rep(NA, 8))
bg_medians <- c(bg_summary[3,1:7], bg_summary[3,9], NA, bg_summary[3,8], NA, bg_summary[3, 10:18], NA, bg_summary[3,19], rep(NA, 8))
harris_means <- c(NA, harris_summary[1, 1:4], NA, harris_summary[1, 5], NA,
                  harris_summary[1,10], NA, harris_summary[1,11], NA, NA, harris_summary[1,7:8], NA, harris_summary[1,9], rep(NA, 3), 
                  harris_summary[1,6], NA, harris_summary[1,12:19])
harris_sds <- c(NA, harris_summary[2, 1:4], NA, harris_summary[2, 5], NA,
                harris_summary[2,10], NA, harris_summary[2,11], NA, NA, harris_summary[2,7:8], NA, harris_summary[2,9], rep(NA, 3), 
                harris_summary[2,6], NA, harris_summary[2,12:19])
harris_medians <- c(NA, harris_summary[3, 1:4], NA, harris_summary[3, 5], NA,
                    harris_summary[3,10], NA, harris_summary[3,11], NA, NA, harris_summary[3,7:8], NA, harris_summary[3,9], rep(NA, 3), 
                    harris_summary[3,6], NA, harris_summary[3,12:19])

sum.table <- cbind(agency_means, county_means, bg_means, harris_means,
                   agency_sds, county_sds, bg_sds, harris_sds,
                   agency_medians, county_medians, bg_medians, harris_medians)

row.names(sum.table) <- c('Crime rate','Murder rate','Robbery rate',
                          'Assault rate','Burglary rate','Larceny rate',
                          'Vehicle theft rate',
                          'Aid value (BG)', 'Aid value (HPBM)', 'Aid quantity (BG)',
                          'Aid quantity (HPBM)', 'Percent poverty', 'Median income', 'Unemployment rate', 'Population', 
                          'Share males','Share blacks','Share age 15–19',
                          'Share age 20–24','Share age 25–34', 'Per cap income', 
                          "BG IV", "HPBM IV 1", "HPBM IV 2", 
                          "HPBM IV 3", "HPBM IV 4", "HPBM IV 5", 
                          "HPBM IV 6", "HPBM IV 7", "HPBM IV 8")
colnames(sum.table) <- rep(c("Agency", "County", "BG", "HPBM"), 3)

sum.table <- xtable(sum.table, digits=1)

#Exporting the table
print(sum.table,
      booktabs = TRUE,
      include.rownames = TRUE,
      file = "./tables/summary-stats.tex",
      caption.placement = "top",
      format.args = list(big.mark = ",", decimal.mark = "."))

#Clear working environment
rm(list=ls())


## Create Supplementary Table 21 (Summary statistics 2010-2012 Subset)
#Load in the data
load('data/BG_milit.RData')
load('data/harris.RData')
load('data/city-aggregate.RData')
load('data/county-aggregate.RData')
harris <- df

#Filter data to 2010-2012 subset
bg_sub <- filter(BG_milit, year > 2009 & year < 2013)
harris_sub <- filter(harris, year > 2009 & year < 2013)
agency_sub <- filter(agency, year > 2009 & year < 2013)
county_sub <- filter(county, year > 2009 & year < 2013)

#Summary stats function
my.summary = function(x, na.rm=TRUE){
  result = c(Mean=mean(x, na.rm=na.rm),
             SD=sd(x, na.rm=na.rm),
             Median=median(x, na.rm=na.rm))
}

#Getting the summary stats for BG
bg_sub <- mutate(bg_sub, shareage2534 = shareage2529+shareage3034)
bg_ind <- c("crime", "murder", "robbery", "assault", "burglary", 
            "larceny", "mvtheft", "tot_q", "total_raw",           
            "PovertyPercentAllAges", "MedianHouseholdIncome", "Unempl_Rate", "pop", "sharemale",            
            "shareblack", "shareage1519", "shareage2024", "shareage2534", "milex_iv")
bg_summary <- as.data.frame(sapply(bg_sub[, bg_ind], my.summary))

#Getting the summary stats for Harris
harris_sub <- mutate(harris_sub, sum_items = exp(log_sum_items))
harris_sub <- mutate(harris_sub, sum_items = if_else(sum_items==1, 0, sum_items))
harris_sub <- mutate(harris_sub, sum_value = exp(log_sum_value_ttl))
harris_sub <- mutate(harris_sub, sum_value = if_else(sum_value==1, 0, sum_value))
harris_ind <- c("Amurders", "Arobbery", "Aassault", "Aburglar", "Avehicle", "PCPI", 
                "UnemploymentRate", "allpop", "blackshare", "sum_value", "sum_items", "zHUdI1",           
                "zHUdI6", "zHUdIland", "zHUdIhidta", "zHUd1", "zHUd6",            
                "zHUdland", "zHUdhitda")
harris_summary <- as.data.frame(sapply(harris_sub[, harris_ind], my.summary))

#Getting the summary stats for Emory Jurisdiction (Agency)
agency_sub <- mutate(agency_sub, shareblack = percblack/100)
agency_ind <- c("allcrime_rate", "Murder_rate", "Robbery_rate", "Assault_rate", "Burglary_rate",           
                "Larceny_rate", "MV_theft_rate", "sumallLESOaidBG", "sumallLESOaidHARRIS",     
                "sumallLESOquantityBG", "sumallLESOquantityHARRIS", "percentinpoverty", "medianhouseholdincome", "unemployment",             
                "Total_Population", "sharemales", "shareblack", "share1519", "share2024",             
                "share2534", "PCPI", "milex_iv", "zHUdI1",           
                "zHUdI6", "zHUdIland", "zHUdIhidta", "zHUd1", "zHUd6",            
                "zHUdland", "zHUdhitda")
agency_summary <- as.data.frame(sapply(agency_sub[, agency_ind], my.summary))

#Getting the summary stats for Emory County
county_sub <- mutate(county_sub, shareblack = percblack/100)
county_ind <- c("allcrime_countysumrate", "MURDER_countysumrate", "ROBBERY_countysumrate", "AGASSLT_countysumrate", "BURGLRY_countysumrate",        
                "LARCENY_countysumrate", "MVTHEFT_countysumrate", "sumallLESOaidBG", "sumallLESOaidHARRIS",      
                "sumallLESOquantityBG", "sumallLESOquantityHARRIS", "percentinpoverty", "medianhouseholdincome", "unemployment",                 
                "Population", "sharemales", "shareblack", "share1519", "share2024",                 
                "share2534", "PCPI", "milex_iv", "zHUdI1",           
                "zHUdI6", "zHUdIland", "zHUdIhidta", "zHUd1", "zHUd6",            
                "zHUdland", "zHUdhitda")
county_summary <- as.data.frame(sapply(county_sub[, county_ind], my.summary))

#Making the table
agency_means <- c(agency_summary[1,])
agency_sds <- c(agency_summary[2,])
agency_medians <- c(agency_summary[3,])
county_means <- c(county_summary[1,])
county_sds <- c(county_summary[2,])
county_medians <- c(county_summary[3,])
bg_means <- c(bg_summary[1,1:7], bg_summary[1,9], NA, bg_summary[1,8], NA, bg_summary[1, 10:18], NA, bg_summary[1, 19], rep(NA, 8))
bg_sds <- c(bg_summary[2,1:7], bg_summary[2,9], NA, bg_summary[2,8], NA, bg_summary[2, 10:18], NA, bg_summary[2,19], rep(NA, 8))
bg_medians <- c(bg_summary[3,1:7], bg_summary[3,9], NA, bg_summary[3,8], NA, bg_summary[3, 10:18], NA, bg_summary[3,19], rep(NA, 8))
harris_means <- c(NA, harris_summary[1, 1:4], NA, harris_summary[1, 5], NA,
                  harris_summary[1,10], NA, harris_summary[1,11], NA, NA, harris_summary[1,7:8], NA, harris_summary[1,9], rep(NA, 3), 
                  harris_summary[1,6], NA, harris_summary[1,12:19])
harris_sds <- c(NA, harris_summary[2, 1:4], NA, harris_summary[2, 5], NA,
                harris_summary[2,10], NA, harris_summary[2,11], NA, NA, harris_summary[2,7:8], NA, harris_summary[2,9], rep(NA, 3), 
                harris_summary[2,6], NA, harris_summary[2,12:19])
harris_medians <- c(NA, harris_summary[3, 1:4], NA, harris_summary[3, 5], NA,
                    harris_summary[3,10], NA, harris_summary[3,11], NA, NA, harris_summary[3,7:8], NA, harris_summary[3,9], rep(NA, 3), 
                    harris_summary[3,6], NA, harris_summary[3,12:19])

sum.table <- cbind(agency_means, county_means, bg_means, harris_means,
                   agency_sds, county_sds, bg_sds, harris_sds,
                   agency_medians, county_medians, bg_medians, harris_medians)

row.names(sum.table) <- c('Crime rate','Murder rate','Robbery rate',
                          'Assault rate','Burglary rate','Larceny rate',
                          'Vehicle theft rate',
                          'Aid value (BG)', 'Aid value (HPBM)', 'Aid quantity (BG)',
                          'Aid quantity (HPBM)', 'Percent poverty', 'Median income', 'Unemployment rate', 'Population', 
                          'Share males','Share blacks','Share age 15–19',
                          'Share age 20–24','Share age 25–34', 'Per cap income', 
                          "BG IV", "HPBM IV 1", "HPBM IV 2", 
                          "HPBM IV 3", "HPBM IV 4", "HPBM IV 5", 
                          "HPBM IV 6", "HPBM IV 7", "HPBM IV 8")
colnames(sum.table) <- rep(c("Agency", "County", "BG", "HPBM"), 3)

sum.table <- xtable(sum.table, digits=1)

#Exporting the table
print(sum.table,
      booktabs = TRUE,
      include.rownames = TRUE,
      file = "./tables/summary-stats-subset.tex",
      caption.placement = "top",
      format.args = list(big.mark = ",", decimal.mark = "."))

#Clear working environment
rm(list=ls())


##F-statistics
library(plm) #for panel data
library(lfe) #for felm function
library(tidyverse)
library(xtable)
library(rlang) #need for functions
source("lib/model-functions-BG.R") #these have all the different models we need to run
source("lib/model-functions-HPBM.R")
year_subset <- FALSE
city_outliers_eliminated <- TRUE

# load data
load('data/city-aggregate.RData')
load('data/county-aggregate.RData')
load('data/BG_milit.RData') # BG data for exact replication
load(file='data/harris.RData') #Harris data for exact replication
harris <- df
rm(df)
## load controlled data
#load('data/city-aggregate-controlled.RData')
#load('data/county-aggregate-controlled.RData')

# restrict data
if(year_subset){
  agency <- agency[agency$year >=2010 & agency$year<=2012,]
  county <- county[county$year >=2010 & county$year <=2012,]
  BG_milit <- BG_milit[BG_milit$year>=2010 & BG_milit$year<=2012,]
  harris <- harris[harris$year>=2010 & harris$year<=2012,]
} 

if(city_outliers_eliminated){
  agency <- agency %>% filter(allcrime_rate < 1000000)
}


# create panel data
agency <- pdata.frame(agency,index=c("g","year"))
county <- pdata.frame(county,index=c("g","year"))
#make it panel and make variables factors
BG_milit$State <- as.factor(BG_milit$State)
BG_milit$fips <- as.factor(BG_milit$fips)

#make data frame panel
BG_milit <- pdata.frame(BG_milit, index = c("fips", "year"))
harris <- pdata.frame(harris, index = c("fips", "year")) 

# run models ------
homicide.items.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'items')
homicide.value.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'value')
homicide.county.items.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'items')
homicide.county.value.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'value')
homicide.items.harris <- run.IV.Harris(Amurders, harris, 'items')
homicide.value.harris <- run.IV.Harris(Amurders, harris, 'value')

allcrime.agency <- run.IV.model.agency(allcrime_rate, agency)
allcrime.county <- run.IV.model.county(allcrime_countysumrate, county)
allcrime.BG <- run.IV.model.BG(crime,BG_milit)
# BG data from stata --------------
#these are from crime estimates
log_sum_aid_BG_est <- c(7.905612,
                        summary(allcrime.agency)$coefficients["`log(1 + sumallLESOaidBG_lag)(fit)`",1],
                        -20.81063,
                        summary(allcrime.county)$coefficients["`log(1 + sumallLESOaidBG_lag)(fit)`",1],
                        -59.29338,
                        summary(allcrime.BG)$coefficients["`ltotal_cost(fit)`",1])

log_sum_aid_BG_se <- c(16.58208,
                       summary(allcrime.agency)$coefficients["`log(1 + sumallLESOaidBG_lag)(fit)`",2],
                       11.99693,
                       summary(allcrime.county)$coefficients["`log(1 + sumallLESOaidBG_lag)(fit)`",2],
                       14.03974,
                       summary(allcrime.BG)$coefficients["`ltotal_cost(fit)`",2])

milex_iv_est <- c(41.28,
                  summary(allcrime.agency$stage1)$coefficients["milex_iv",1],
                  46.64,
                  summary(allcrime.county$stage1)$coefficients["milex_iv",1],
                  17.59,
                  summary(allcrime.BG$stage1)$coefficients["milex_iv",1])

milex_iv_se <- c(4.511,
                 summary(allcrime.agency$stage1)$coefficients["milex_iv",2],
                 3.472,
                 summary(allcrime.county$stage1)$coefficients["milex_iv",2],
                 2.510,
                 summary(allcrime.BG$stage1)$coefficients["milex_iv",2])


n_obs <- c(44963,
           summary(allcrime.agency)$N,
           15683,
           summary(allcrime.county)$N,
           17807,
           summary(allcrime.BG)$N)

group <- c("Agency","Agency", 
           "County_all_crime", "County_all_crime",
           "BG_all_crime","BG_all_crime")

software <- c("Stata","R","Stata","R","Stata","R")

KP_wald_f <- c(83.760,NA,180.500,NA, 49.113,NA) #fixing clustering to state level fixed f test
CD_wald_f <- c(1066.044,NA, 749.330,NA,  258.545,NA)
cluster_wald_f <- c(NA,
                    waldtest(allcrime.agency$stage1, R = 'milex_iv', type = 'cluster')["F"],
                    NA,
                    waldtest(allcrime.county$stage1, R = 'milex_iv', type = 'cluster')["F"],
                    NA,
                    waldtest(allcrime.BG$stage1, R = 'milex_iv', type = 'cluster')["F"])


f_stats <- data.frame(log_sum_aid_BG_est,
                      log_sum_aid_BG_se,
                      milex_iv_est,
                      milex_iv_se,
                      n_obs,
                      group,
                      software,
                      KP_wald_f,
                      CD_wald_f,
                      cluster_wald_f)

f_stats

tab1 <- xtable(f_stats)
print(tab1, type = 'latex',file='tables/f_stats-BG.tex')

## Now do HPBM ---------
## wald test in R ----
HPBM_IVs_item <- c("zHUdI1_lag", "zHUdIland_lag", "zHUdI6_lag","zHUdIhidta_lag")
w1 <- waldtest(homicide.items.agency.Harris$stage1, R = HPBM_IVs_item, type = 'cluster')["F"]
w2 <- waldtest(homicide.county.items.harris$stage1, R = HPBM_IVs_item, type = 'cluster')["F"]
w3 <- waldtest(homicide.items.harris$stage1, R = HPBM_IVs_item, type = 'cluster')["F"]
# data from stata
#these are from homicide estimates items
log_sum_items_HPBM_est <- c(-.8756096,
                            summary(homicide.items.agency.Harris)$coefficients["`log(1 + sumallLESOquantityHARRIS_lag)(fit)`",1],
                            -.2834097,
                            summary(homicide.county.items.harris)$coefficients["`log(1 + sumallLESOquantityHARRIS_lag)(fit)`",1],
                            -.2206904,
                            summary(homicide.items.harris)$coefficients["`log_sum_items_lag(fit)`",1])

log_sum_items_HPBM_se <- c(1.127337,
                           summary(homicide.items.agency.Harris)$coefficients["`log(1 + sumallLESOquantityHARRIS_lag)(fit)`",2],
                           .5093377,
                           summary(homicide.county.items.harris)$coefficients["`log(1 + sumallLESOquantityHARRIS_lag)(fit)`",2],
                           .3024331,
                           summary(homicide.items.harris)$coefficients["`log_sum_items_lag(fit)`",2])

iv_est <- c(0.00000254, 0.00000225, 0.00177,-0.000000773,
            summary(homicide.items.agency.Harris$stage1)$coefficients[HPBM_IVs_item,1],
            0.0000290,  0.000000953, -0.00251,0.00000764,
            summary(homicide.county.items.harris$stage1)$coefficients[HPBM_IVs_item,1],
            2.127,8.13e-08,-37.52,0.00000156,
            summary(homicide.items.harris$stage1)$coefficients[HPBM_IVs_item,1])

iv_se <- c(0.00000628, 0.000000208, 0.00113,0.000000447,
           summary(homicide.items.agency.Harris$stage1)$coefficients[HPBM_IVs_item,2],
           0.0000519, 0.000000589, 0.00280, 0.00000148,
           summary(homicide.county.items.harris$stage1)$coefficients[HPBM_IVs_item,2],
           0.708,0.000000102,17.07,0.000000224,
           summary(homicide.items.harris$stage1)$coefficients[HPBM_IVs_item,2])

iv_type <- c(rep(x=HPBM_IVs_item,6))
#rep(c("zHUdI1_lag", "zHUdIland_lag", "zHUdI6_lag","zHUdIhidta_lag"),2)

n_obs <- c( 44693,
            summary(homicide.items.agency.Harris)$N,
            15460,
            summary(homicide.county.items.harris)$N,
            36671,
            summary(homicide.items.harris)$N)

group <- c(rep("Agency",8),
           rep("County",8),
           rep("HPBM",8))
software <- rep(c(rep("Stata",4),rep("R",4)),3)

type = rep("items",6)


CD_wald_f <- c( 41.455,
                NA,
                13.599,
                NA,
                149.212,
                NA)
KP_wald_f <- c(25.608,
               NA,
               8.512,
               NA,
               16.132,
               NA)
cluster_wald_f <- c(NA,
                    w1,
                    NA,
                    w2,
                    NA,
                    w3)


iv <- data.frame(iv_est, iv_se, iv_type,group,software)
test <- pivot_wider(iv,names_from = iv_type, values_from=c(iv_est,iv_se))

f_stats <- data.frame(log_sum_items_HPBM_est,
                      log_sum_items_HPBM_se,
                      test,
                      n_obs,
                      KP_wald_f,
                      CD_wald_f,
                      cluster_wald_f)

View(f_stats)

tab1 <- xtable(f_stats)
print(tab1, type = 'latex',file='tables/f_stats-HPBM-item.tex')



## ---------------------------------------------------- ##
## Repeat but for value of items ----
## ---------------------------------------------------- ##
HPBM_IVs_value <- c("zHUd1_lag", "zHUd6_lag", "zHUdland_lag", "zHUdhitda_lag")
wv1 <- waldtest(homicide.value.agency.Harris$stage1, R = HPBM_IVs_value, type = 'cluster')["F"]
wv2 <- waldtest(homicide.county.value.harris$stage1, R = HPBM_IVs_value, type = 'cluster')["F"]
wv3 <- waldtest(homicide.value.harris$stage1, R = HPBM_IVs_value, type = 'cluster')["F"]
# data from stata
#these are from homicide estimates value
log_sum_value_HPBM_est <- c( .1446227,
                             summary(homicide.value.agency.Harris)$coefficients["`log(1 + sumallLESOaidHARRIS_lag)(fit)`",1],
                             .1374266,
                             summary(homicide.county.value.harris)$coefficients["`log(1 + sumallLESOaidHARRIS_lag)(fit)`",1],
                             .3487729,
                             summary(homicide.value.harris)$coefficients["`log_sum_value_ttl_lag(fit)`",1])

log_sum_value_HPBM_se <- c(.2959486,
                           summary(homicide.value.agency.Harris)$coefficients["`log(1 + sumallLESOaidHARRIS_lag)(fit)`",2],
                           .1782375,
                           summary(homicide.county.value.harris)$coefficients["`log(1 + sumallLESOaidHARRIS_lag)(fit)`",2],
                           .2818826,
                           summary(homicide.value.harris)$coefficients["`log_sum_value_ttl_lag(fit)`",2])

iv_est <- c(  4.31e-10,  -0.000000333,1.18e-09, 9.45e-10,
              summary(homicide.value.agency.Harris$stage1)$coefficients[HPBM_IVs_value,1],
              4.69e-09, -0.00000879, -5.29e-11, 6.46e-09,
              summary(homicide.county.value.harris$stage1)$coefficients[HPBM_IVs_value,1],
              9.058, -203.8, 0.0790, 0.138,
              summary(homicide.value.harris$stage1)$coefficients[HPBM_IVs_value,1])

iv_se <- c(2.41e-09,0.000000703,1.12e-10,2.68e-10,
           summary(homicide.value.agency.Harris$stage1)$coefficients[HPBM_IVs_value,2],
           6.41e-08,0.00000301,6.03e-10,1.46e-09,
           summary(homicide.county.value.harris$stage1)$coefficients[HPBM_IVs_value,2],
           2.715,120.2,0.0671,0.138,
           summary(homicide.value.harris$stage1)$coefficients[HPBM_IVs_value,2])

iv_type <- rep(HPBM_IVs_value,6)

n_obs <- c(44693,
           summary(homicide.value.agency.Harris)$N,
           15460,
           summary(homicide.county.value.harris)$N,
           36671,
           summary(homicide.value.harris)$N)

group <- c(rep("Agency",8),
           rep("County",8),
           rep("HPBM",8))

software <- rep(c(rep("Stata",4),rep("R",4)),3)
CD_wald_f <- c( 50.880,
                NA,
                11.226,
                NA,
                20.528,
                NA)
KP_wald_f <- c( 28.862,
                NA,
                7.750,
                NA,
                4.145,
                NA)

cluster_wald_f <- c(NA,
                    wv1,
                    NA,
                    wv2,
                    NA,
                    wv3)

iv <- data.frame(iv_est, iv_se, iv_type,group,software)
test <- pivot_wider(iv,names_from = iv_type, values_from=c(iv_est,iv_se))

f_stats <- data.frame(log_sum_value_HPBM_est,
                      log_sum_value_HPBM_se,
                      test,
                      n_obs,
                      KP_wald_f,
                      CD_wald_f,
                      cluster_wald_f)
View(f_stats)


tab2 <- xtable(f_stats)
print(tab2, type = 'latex',file='tables/f_stats-HPBM-value.tex')

## ------------------------------------------------------- ##
## Repeat on subset of data fromo 2010-2012 for overlap time periods ----
## ------------------------------------------------------- ##
source("lib/model-functions-BG.R") #these have all the different models we need to run
source("lib/model-functions-HPBM.R")
year_subset <- TRUE
city_outliers_eliminated <- TRUE

# load data
load('data/city-aggregate.RData')
load('data/county-aggregate.RData')
load('data/BG_milit.RData') # BG data for exact replication
load(file='data/harris.RData') #Harris data for exact replication
harris <- df
rm(df)
## load controlled data
#load('data/city-aggregate-controlled.RData')
#load('data/county-aggregate-controlled.RData')

# restrict data
if(year_subset){
  agency <- agency[agency$year >=2010 & agency$year<=2012,]
  county <- county[county$year >=2010 & county$year <=2012,]
  BG_milit <- BG_milit[BG_milit$year>=2010 & BG_milit$year<=2012,]
  harris <- harris[harris$year>=2010 & harris$year<=2012,]
} 

if(city_outliers_eliminated){
  agency <- agency %>% filter(allcrime_rate < 1000000)
}

# create panel data
agency <- pdata.frame(agency,index=c("g","year"))
county <- pdata.frame(county,index=c("g","year"))

#make it panel and make variables factors
BG_milit$State <- as.factor(BG_milit$State)
BG_milit$fips <- as.factor(BG_milit$fips)

#make data frame panel
BG_milit <- pdata.frame(BG_milit, index = c("fips", "year"))
harris <- pdata.frame(harris, index = c("fips", "year")) 

# run models ------
homicide.items.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'items')
homicide.value.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'value')
homicide.county.items.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'items')
homicide.county.value.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'value')
homicide.items.harris <- run.IV.Harris(Amurders, harris, 'items')
homicide.value.harris <- run.IV.Harris(Amurders, harris, 'value')

allcrime.agency <- run.IV.model.agency(allcrime_rate, agency)
allcrime.county <- run.IV.model.county(allcrime_countysumrate, county)
allcrime.BG <- run.IV.model.BG(crime,BG_milit)

# BG subset wald test ----
w_BG_agency_subset <- waldtest(allcrime.agency$stage1, R = 'milex_iv', type = 'cluster')["F"]
w_BG_county_subset <- waldtest(allcrime.county$stage1, R = 'milex_iv', type = 'cluster')["F"]
w_BG_subset <- waldtest(allcrime.BG$stage1, R = 'milex_iv', type = 'cluster')["F"]

waldBG <- data.frame(BG_agency = w_BG_agency_subset,
                     BG_county = w_BG_county_subset,
                     BG = w_BG_subset)

xtable(waldBG)
# HPBM subset wald test ----
HPBM_IVs_item <- c("zHUdI1_lag", "zHUdIland_lag", "zHUdI6_lag","zHUdIhidta_lag")
w_HPBM_agency_subset_item <- waldtest(homicide.items.agency.Harris$stage1, R = HPBM_IVs_item, type = 'cluster')["F"]
w_HPBM_county_subset_item <- waldtest(homicide.county.items.harris$stage1, R = HPBM_IVs_item, type = 'cluster')["F"]
w_HPBM_subset_item <- waldtest(homicide.items.harris$stage1, R = HPBM_IVs_item, type = 'cluster')["F"]
waldHPBM_item <- data.frame(HPBM_agency = w_HPBM_agency_subset_item,
                            HPBM_county = w_HPBM_county_subset_item,
                            HPBM = w_HPBM_subset_item)

HPBM_IVs_value <- c("zHUd1_lag", "zHUd6_lag", "zHUdland_lag", "zHUdhitda_lag")
w_HPBM_agency_subset <- waldtest(homicide.value.agency.Harris$stage1, R = HPBM_IVs_value, type = 'cluster')["F"]
w_HPBM_county_subset <- waldtest(homicide.county.value.harris$stage1, R = HPBM_IVs_value, type = 'cluster')["F"]
w_HPBM_subset <- waldtest(homicide.value.harris$stage1, R = HPBM_IVs_value, type = 'cluster')["F"]
waldHPBM <- data.frame(HPBM_agency = w_HPBM_agency_subset,
                       HPBM_county = w_HPBM_county_subset,
                       HPBM = w_HPBM_subset)
xtable(waldHPBM)

rm(list=ls())

